function [ll, llgrad] = Phat(theta,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow)

N=size(X,1);
KZ=size(Z,2)/5;
NS=size(V_R,2);
nBeta=length(reb{1});
nPhi=length(mpow{1});
nGamma=length(reb{2});
nChi=length(mpow{2});
nYY=size(YY,1);

Xrebel=X(:,reb{1}); Xmpow=X(:,mpow{1});

ZUS=Z(:,1:KZ);
ZUK=Z(:,KZ+(1:KZ)); 
ZFRN=Z(:,2*KZ+(1:KZ));
ZRUS=Z(:,3*KZ+(1:KZ));
ZCHN=Z(:,4*KZ+(1:KZ));

ZrebelUS=ZUS(:,reb{2}); ZmpowUS=ZUS(:,mpow{2});
ZrebelUK=ZUK(:,reb{2}); ZmpowUK=ZUK(:,mpow{2});
ZrebelFRN=ZFRN(:,reb{2}); ZmpowFRN=ZFRN(:,mpow{2});
ZrebelRUS=ZRUS(:,reb{2}); ZmpowRUS=ZRUS(:,mpow{2});
ZrebelCHN=ZCHN(:,reb{2}); ZmpowCHN=ZCHN(:,mpow{2});

beta=theta(1:nBeta,1);

phi_US=theta(nBeta+[1,5+(1:nPhi-1)],1);
phi_UK=theta(nBeta+[2,5+(1:nPhi-1)],1);
phi_France=theta(nBeta+[3,5+(1:nPhi-1)],1);
phi_Russia=theta(nBeta+[4,5+(1:nPhi-1)],1);
phi_China=theta(nBeta+[5,5+(1:nPhi-1)],1);
    
gamma_US=theta(nBeta+4+nPhi+1,1);
gamma_UK=theta(nBeta+4+nPhi+2,1);
gamma_France=theta(nBeta+4+nPhi+3,1);
gamma_Russia=theta(nBeta+4+nPhi+4,1);
gamma_China=theta(nBeta+4+nPhi+5,1);
gamma0=theta(nBeta+4+nPhi+5+(1:nGamma),1);

chi=theta(nBeta+4+nPhi+5+nGamma+(1:nChi),1);

delta_US_UK=theta(nBeta+4+nPhi+5+nGamma+nChi+1,1);
delta_US_France=theta(nBeta+4+nPhi+5+nGamma+nChi+2,1);
delta_US_Russia=theta(nBeta+4+nPhi+5+nGamma+nChi+3,1);
delta_US_China=theta(nBeta+4+nPhi+5+nGamma+nChi+4,1);
delta_UK_France=theta(nBeta+4+nPhi+5+nGamma+nChi+5,1);
delta_UK_Russia=theta(nBeta+4+nPhi+5+nGamma+nChi+6,1);
delta_UK_China=theta(nBeta+4+nPhi+5+nGamma+nChi+7,1);
delta_France_Russia=theta(nBeta+4+nPhi+5+nGamma+nChi+8,1);
delta_France_China=theta(nBeta+4+nPhi+5+nGamma+nChi+9,1);
delta_Russia_China=theta(nBeta+4+nPhi+5+nGamma+nChi+10,1);

lambda=theta(nBeta+4+nPhi+5+nGamma+nChi+11:end,1);

ll=0;
llgrad=0;

for n=1:N
    
    u_R = YY(:,1) .* (Xrebel(n,:)*beta+YY(:,2:end)*[gamma_US+ZrebelUS(n,:)*gamma0;gamma_UK+ZrebelUK(n,:)*gamma0;...
        gamma_France+ZrebelFRN(n,:)*gamma0;gamma_Russia+ZrebelRUS(n,:)*gamma0;gamma_China+ZrebelCHN(n,:)*gamma0]);
    
    u_US = YY(:,2) .* (Xmpow(n,:)*phi_US + ZmpowUS(n,:)*chi +...
        YY(:,3:6) * [delta_US_UK; delta_US_France; delta_US_Russia; delta_US_China]);
    
    u_UK = YY(:,3) .* (Xmpow(n,:)*phi_UK + ZmpowUK(n,:)*chi +...
        YY(:,[2,4:6])*[delta_US_UK; delta_UK_France; delta_UK_Russia; delta_UK_China]);
    
    u_France = YY(:,4) .* (Xmpow(n,:)*phi_France + ZmpowFRN(n,:)*chi +...
        YY(:,[2:3,5:6])*[delta_US_France; delta_UK_France; delta_France_Russia; delta_France_China]);
    
    u_Russia = YY(:,5) .* (Xmpow(n,:)*phi_Russia + ZmpowRUS(n,:)*chi +...
        YY(:,[2:4,6])*[delta_US_Russia; delta_UK_Russia; delta_France_Russia; delta_Russia_China]);
    
    u_China = YY(:,6) .* (Xmpow(n,:)*phi_China + ZmpowCHN(n,:)*chi +...
        YY(:,2:5)*[delta_US_China; delta_UK_China; delta_France_China; delta_Russia_China]);
       
    V=[V_R(:,:,n);V_US(:,:,n);V_UK(:,:,n);V_France(:,:,n);V_Russia(:,:,n);V_China(:,:,n)]';
    U=[u_R;u_US;u_UK;u_France;u_Russia;u_China]';
    U0n=reshape(U0(:,:,n),6*nYY,1)';
    weights = exp(-sum((V-U).^2,2)/2 + sum((V-U0n).^2,2)/2) / NS;
    
    
    if nargout<2
        
        Pra = zeros(NS, 1);
        for b=1:NS
            Pra(b,1) = profprob(Y(n,1), DMES{b,n}, intervEq{b,n}, V_R(1,b,n), lambda);
        end
    
        ll = ll - log(weights'*Pra);

    else
        
        %Gradients w.r.t. payoffs
        Du=zeros(6*nYY,length(theta)-length(lambda));
        % rebel payoffs  
        Du(1:nYY,[1:nBeta,nBeta+4+nPhi+(1:5+nGamma)])=[YY(:,1)*Xrebel(n,:),YY(:,2:end),YY(:,2)*ZrebelUS(n,:)+...
            YY(:,3)*ZrebelUK(n,:)+YY(:,4)*ZrebelFRN(n,:)+YY(:,5)*ZrebelRUS(n,:)+YY(:,6)*ZrebelCHN(n,:)];
        % US payoffs
        Du(nYY+(1:nYY),[nBeta+[1,5+(1:nPhi-1)],nBeta+4+nPhi+5+nGamma+(1:nChi),nBeta+4+nPhi+5+nGamma+nChi+(1:4)])=YY(:,2).*[YY(:,1)*[Xmpow(n,:),ZmpowUS(n,:)],YY(:,3:6)];
        % UK
        Du(2*nYY+(1:nYY),[nBeta+[2,5+(1:nPhi-1)],nBeta+4+nPhi+5+nGamma+(1:nChi),nBeta+4+nPhi+5+nGamma+nChi+[1,5:7]])=YY(:,3).*[YY(:,1)*[Xmpow(n,:),ZmpowUK(n,:)],YY(:,[2,4:6])];
        % France 
        Du(3*nYY+(1:nYY),[nBeta+[3,5+(1:nPhi-1)],nBeta+4+nPhi+5+nGamma+(1:nChi),nBeta+4+nPhi+5+nGamma+nChi+[2,5,8:9]])=YY(:,4).*[YY(:,1)*[Xmpow(n,:),ZmpowFRN(n,:)],YY(:,[2:3,5:6])];
        % Russia
        Du(4*nYY+(1:nYY),[nBeta+[4,5+(1:nPhi-1)],nBeta+4+nPhi+5+nGamma+(1:nChi),nBeta+4+nPhi+5+nGamma+nChi+[3,6,8,10]])=YY(:,5).*[YY(:,1)*[Xmpow(n,:),ZmpowRUS(n,:)],YY(:,[2:4,6])];
        % China
        Du(5*nYY+(1:nYY),[nBeta+[5,5+(1:nPhi-1)],nBeta+4+nPhi+5+nGamma+(1:nChi),nBeta+4+nPhi+5+nGamma+nChi+[4,7,9:10]])=YY(:,6).*[YY(:,1)*[Xmpow(n,:),ZmpowCHN(n,:)],YY(:,2:5)];

        Pra = zeros(NS, 1);
        Grad = zeros(NS, length(lambda));
        for b=1:NS
            [Pra(b,1), Grad(b,:)] = profprob(Y(n,1), DMES{b,n}, intervEq{b,n}, V_R(1,b,n), lambda);
        end
    
        ln = weights'*Pra;
        ll = ll - log(ln);
        llgrad = llgrad - weights'*[Pra.*((V-U)*Du),Grad]/ln;     
        
    end         
    
end


